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Multi-Level  Approximation  to  Scattered 
Data  Using  Inverse  Multiquadrics 

S.  J.  Hales  and  J.  Levesley 


Abstract.  A method  of  finding  local  approximations  is  used  to  thin 
data  before  a hierarchical  iterative  refinement  scheme  is  employed  in  con- 
junction with  domain  decomposition.  The  interpolation  problem  on  each 
sub-domain  is  solved  by  using  the  same  stored  inverse.  The  approxima- 
tion power  of  the  inverse  multiquadric  is  exploited  whilst  overcoming  the 
computational  difficulties  associated  with  globally  supported  basis  func- 
tions. 


§1.  Introduction 

Radial  basis  functions  have  been  widely  used  for  multivariate  interpolation  of 
scattered  data,  see  [4]  for  a summary.  An  interpolant  is  generated  by  a linear 
combination  of  basis  functions  cj>  at  distinct  centres  x*,  i = 1 

N 

«(*)  = d) 

i=l 

constrained  by  s(xj)  = /,  , * = 1, ...,1V,  where  T : 1R“  >—>  ]R  and  = 
T(xi).  The  interpolation  matrix  A € IRiV  x 1R‘V  is  given  by  Ai<jj<jv  = 
0(||  a Zi  — ac j- 1|),  and  A satisfies 

AA  = /,  (2) 

where  A = [Ai  • • ■ \n}t  and  / = [/i  • • • /jv]T- 

Common  choices  for  <j>  in  this  setting  are  given  in  [6], 

0(||x  — a^i  ||)  = exp(— c2||x  — Xi||2),  Gaussian, 

<f)(\\x  — Xj||)  = (c2  + ||x  — XiH2)-1^2,  Inverse  multiquadric, 

where  c is  a constant  shape  parameter.  With  a small  modification  to  the 
scheme,  the  thin  plate  spline  and  multiquadric  are  also  used. 
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The  above  parameter-dependent  functions  are  good  at  approximating 
data  for  certain  values  of  c,  but  these  cause  inherent  ill-conditioning  in  A. 
Schaback  [7]  explains  this  phenomenon  by  means  of  an  “Uncertainty  Relation” 
between  upper  bounds  on  errors  for  interpolants  of  the  form  (1),  and  lower 
bounds  on  the  smallest  eigenvalue  of  A.  Iterative  techniques  for  solving  such 
badly  conditioned  systems  often  suffer  from  poor  rates  of  convergence,  and 
therefore  computationally  expensive  direct  methods  have  to  be  employed. 

In  the  inverse  multiquadric  case,  large  values  of  c achieve  good  initial  ap- 
proximations to  smooth  data,  whilst  smaller  values  produce  functions  capable 
of  resolving  fine  detail.  Ideally,  such  properties  could  be  exploited  without 
having  to  solve  (2)  directly. 

Since  s is  evaluated  at  points  y ^ aq  where  an  approximation  is  required, 
a global  solution  incorporating  all  N centres  may  be  inappropriate.  Further, 
it  is  unnecessary  to  find  A such  that  jj/  — AA]^  <C  |.F(£/)  — s(j/)|,  since  the 
accuracy  of  s(y ) is  limited  by  the  approximating  power  of  <fi.  Rather  than 
searching  for  a complete  global  solution,  this  suggests  that  attention  may  be 
focussed  on  small  regions  around  evaluation  points.  Moreover,  the  aim  is 
to  obtain  a solution  such  that  the  residual  and  approximation  accuracy  are 
comparable,  for  little  is  to  be  gained  by  having  a small  residual,  while  the 
approximation  power  of  the  basis  functions  limits  the  final  accuracy. 

In  Section  2,  local  approximations  are  used  to  convert  irregular  data  to 
a regular  mesh  of  approximate  function  values.  Whilst  the  method  can  be 
generalised  to  Rd,  the  description  and  examples  are  given  in  1R2.  The  system 
of  equations  associated  with  the  gridded  data  is  inverted  and  used  to  solve 
subsequent  systems. 

Floater  & Iske  [1]  demonstrate  the  benefits  of  a multi-level  approach 
to  approximation,  and  the  theoretical  foundation  is  provided  by  Narcowich, 
Schaback  & Ward  [5].  Section  3 describes  the  present  hierarchical  iterative 
refinement  algorithm,  and  explains  the  computational  advantages  of  domain 
decomposition  and  the  use  of  a stored  inverse. 


§2.  Local  Solutions  and  Gridding  Data 

If  the  function  T is  not  arbitrary,  but  arises  from  a physical  system,  then  some 
degree  of  smoothness  can  be  assumed.  A smooth  data  set  can  be  significantly 
thinned  whilst  retaining  general  information  about  its  behaviour.  Floater  & 
Iske  [1]  demonstrate  that  Delaunay  triangulation  can  be  used  to  optimise  the 
uniformity  of  data,  and  provide  a good  thinning  algorithm.  Such  triangulation 
and  assembling  of  data  is  computationally  expensive  for  excessively  large  N. 
An  O(N)  method  of  finding  uniform  approximate  data  is  presented. 

An  approximation  to  T at  a point  y € lRd  is  achieved  by  solving  a small 
interpolation  problem  centred  on  y.  The  closest  q points  in  X to  y are  inter- 
polated by  inverse  multiquadrics  with  shape  parameter  C(oca;,  and  evaluated 
at  y.  Since  q can  be  as  low  as  20  ~ 30,  c/oca;  can  be  relatively  large  be- 
fore the  matrix  ill-conditioning  becomes  unacceptable,  thus  yielding  a good 
approximation. 
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This  method  is  highly  parallelizable,  and  large  data-sets  can  be  dealt 
with  without  the  need  of  assembling  or  storing  the  matrix  A. 

Finding  the  optimum  shape  parameter  on  regular  or  scattered  data  re- 
mains an  open  problem,  as  shown  in  [3].  There  is  no  obvious  correlation 
between  point  spacing  and  a good  choice  of  c;oco;.  The  best  shape  parameter 
is  generally  found  by  increasing  the  value  of  cjocoj  until  just  prior  to  machine 
precision  breakdown. 

Let  Y = {yi, . . . , yn  2}  be  the  set  of  points  on  the  nxn  regular  unit  grid. 
If  the  previous  local  approximation  technique  is  applied  to  each  yi,  then  the 
irregular  data  can  be  transformed  to  a regular  grid  with  approximate  function 
values  /j.  The  aim  is  to  find  a global  approximation  using  the  new  data  at 
the  grid  points. 

After  converting  scattered  data  to  a regular  grid,  certain  approximation 
techniques  become  available  which  would  otherwise  have  been  difficult  to  im- 
plement. Polynomial  tensor  product  splines  can  be  efficiently  employed  to 
approximate  a solution  from  the  given  gridded  data.  To  find  such  an  approxi- 
mation at  a point  z , z must  lie  inside  a (d  + 1)  x (d  + 1)  subgrid  of  the  regular 
points,  where  d is  the  degree  of  the  Lagrange  polynomials  to  be  used.  Let 
the  points  of  such  a subgrid  be  labelled  and  have  function  values  fij  for 
i,j  = 1, . . . ,d  + 1.  The  univariate  Lagrange  polynomials  Li(x)  and  iJ (y)  are 
constructed  such  that 


Li(tk.)  = 6 ? and  V(U)  = 6$. 

The  polynomial  tensor  product  spline  faj  is  defined  to  be 

<t>ij{z)  = Ll(z)  ■ Lj(z). 

The  approximation  at  z is  given  by 

d+1  d+ 1 

fa 

i=l  j= 1 

Alternatively,  a thinned  global  interpolant  of  the  form  (1)  can  be  achieved  by 
solving 


Bp  = f, 


(3) 


where  = fa\\yi  — 2/j  1 1 ) > P = \pi"  ' Pn]T  and  / = [/1  • • • f„]T- 

This  amounts  to  finding  an  interpolant  s to  a thinned  approximation  of 
the  initial  data.  The  local  approximation  errors  | f(yi)  — f(Vi) I limit  the  final 
accuracy  of  s. 

The  inverse  of  B need  only  be  computed  once,  and  then  stored  for  future 
use.  All  scattered  data  problems  can  then  be  scaled  and  transformed  to  the 
regular  grid  Y,  whereupon  p is  given  by  the  matrix-vector  product  p = B~xf. 
Only  half  of  the  entries  of  B^1  need  to  be  stored  since  B~l  = (B~1)T . 
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Fig.  1.  An  example  of  approximation  vs.  residual. 

Too  large  a matrix  B causes  storage  problems,  and  difficulties  in  calculat- 
ing the  inverse.  As  c increases,  the  approximation  improves,  but  the  residual 
11/  — B/i ||oo  grows.  A value  for  c is  chosen  before  the  approximation  begins  to 
deteriorate  due  to  the  rise  in  the  residual.  As  an  example,  the  function  T = 1 
is  approximated  on  the  unit  square  using  inverse  multiquadrics  by  interpolat- 
ing fc  = T{yi),  i = 1, . . . ,400  using  (3).  The  approximation  is  evaluated  at 
1000  random  points.  The  results  in  Figure  1 are  typical  for  smooth  functions, 
but  the  consequent  choice  of  c is  only  a guide,  and  does  not  guarantee  success 
for  all  T. 


§3.  Hierarchical  Iterative  Refinement 

The  hierarchical  method  uses  increasingly  dense  subsets  of  X to  refine  the 
current  approximation;  see  [5].  Let  \k  = {iq, . . . ,&Nk}  £ A,  such  that 
Nk+i  > Nk.  Let  St  be  the  current  approximation,  and  r k be  the  full  global 
residual  at  the  kth  level, 


Tk(Xi)  - f{Xi)  - Sk(Xi). 


Let  rk  be  the  kth  residual  over  the  points  in  Xk ■ This  thinned  global 
residual  is  interpolated  by 


Nk 

tk{x)  = ^7i^(||x-ii||),  (4) 

t = l 

where  <t>k{\\x-Xi\\)  = (c£  + ||a;-i1||2)-1''2  , and  tk(xi)  = rk(xt)  , i = 1, . . . , Nk. 
The  initial  interpolant  si  = 0 is  updated  by 


sk+ 1 = sk  + tk. 


(5) 
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The  technique  of  gridding  data  in  Section  2 is  used  to  find  an  approximate 
function  value  for  every  point  in  Yp.  Therefore,  (5)  is  replaced  by  Sfc+i  = 
Sk  + tk,  where  Sk  is  the  current  approximation  to  regular  approximate  data, 

si  = 0. 

The  value  ci  can  be  relatively  large  to  give  a good  initial  approximation. 
As  Nk  increases,  Ck  has  to  be  reduced  to  ensure  computational  solvability.  The 
decrease  in  Ck  introduces  tighter  basis  functions  which  improve  the  resolution 
of  the  approximation. 

A method  of  data  thinning  is  required  to  determine  the  points  in  \k  ■ The 
dense  systems  arising  from  (4)  have  to  be  solved  directly,  but  this  is  imprac- 
tical for  large  Nk-  To  overcome  such  complications,  domain  decomposition  is 
applied  to  each  Xk- 

The  levels  of  the  hierarchy  have  to  be  computed  sequentially,  but  by 
using  domain  decomposition  each  sub-domain  can  be  dealt  with  in  parallel. 
Moreover,  each  such  solution  only  requires  a single  matrix-vector  product. 

To  put  this  in  the  current  context,  each  Xk  is  constructed  from  overlapping 
square  grids  Yp  , where  p = 1, . . . , m*..  These  square  grids  need  not  be  the 
same  size  or  of  similar  orientation,  but  must  contain  an  equal  number  of  points. 
Each  sub-domain  Yp  consists  of  an  inner  region,  where  the  approximation  is 
finally  evaluated,  and  an  overlap.  Special  attention  has  to  be  given  to  sub- 
domains  whose  edges  coincide  with  the  boundary  of  X. 

At  the  kth  level,  m*,  sub-domain  interpolation  problems  need  to  be  solved. 
Since  B is  invariant  under  shifts  and  rotations  of  the  centres  yi,  the  stored 
B~l  can  be  invoked.  If  the  centres  are  scaled  yi  oj/j,  this  amounts  to  a 
change  in  the  shape  parameter. 

Recall  that  Bi<;,j<n  = <i>(\\yi  - S/j  H),  where  y{  € [0, 1]  X [0, 1].  Now, 


<K\\v-Vi\\)  = {<?  + h - Vi\\2)  1/2 

= a(a2c2  + || ay  - ayi\\2)~1/2. 

Let  Wi  = ayi  and  define  i>(\\w  — u>i||)  = a(a2 c2  + ||w  — Wj||2)-1/2.  Then 
Bi <ij<„  = <t>(\\yi~yj\\)  = il){\\wi-Wj\\)  where  w*  € [0,  a]  x [0,  a).  Therefore 
by  using  the  matrix  B,  a new  inverse  multiquadric  is  created  at  scaled  points 
with  shape  parameter  ac. 

Each  of  the  thinned  global  interpolation  problems  (4)  can  be  decomposed 
and  solved  by  multiple  applications  of  the  stored  inverse  B~1.  Continued  use 
of  the  same  inverse  naturally  introduces  tighter  basis  functions  suitable  for 
approximating  typical  residuals. 

§4.  Numerical  Results 

We  give  an  example  where  the  above  scheme  is  used  to  approximate  Franke’s 
function  [2]  over  10000  scattered  points  in  the  unit  square  in  H2: 

X{u,  v ) = o.75e-0-25(9u_2)2-°'25(9^2)2  + o.75e-(9“-2)2/49-(9v-2)2/10 

+ o.5e~°-25(9u-7)2~°-25(9,,-3)2  - 0.2e-(9u-4)2~(9t,~7)2- 
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Level 

No.  of 
domains 
mk 

Shape 

parameter 

Ck 

Overlap 

Max.  error  in 
gridded  data 
Wfi-Moo 

Max.  error  in 
solution 
|j*^*  &k  ||oo 

1 

1 

0.25 

0 

8.241  x 10~5 

1.371  x 10~4 

2 

4 

0.138 

1/36 

2.424  x 10“6 

1.838  x HT5 

3 

16 

0.0688 

1/72 

3.189  x 10~6 

3.189  x 10~6 

Tab.  1.  Error  for  Franke’s  function. 

The  localised  interpolation  problems  are  solved  directly  using  Gauss  Elim- 
ination, with  q = 20  and  C[ocai  = 0.2.  It  is  the  error  function  at  each  level 
which  is  approximated  locally,  and  not  the  original  function  T . The  square 
sub-domain  grids  Yp  are  comprised  of  21  x 21  equally  spaced  points.  For  ease 
of  implementation,  the  sub-domains  used  for  a particular  level  are  of  equal 
size.  The  overlaps  between  sub-domains  therefore  consist  of  one  or  two  mesh 
points,  depending  on  position.  The  key  interpolation  matrix  B is  constructed 
from  inverse  multiquadrics  with  j/j  € [0,  l]2,  c = 0.25,  and  B~l  is  generated  us- 
ing Matlab.  The  domain  decomposition  is  straightforward  on  the  unit  square 
with  m*.  = 4k~1.  The  thinned  global  interpolants  Sk  are  evaluated  at  points 
U € [0,  l]2.  Table  1 shows  the  error  in  the  approximated  data  at  the  regular 
grid  points  \f{yi)  — f(yi) |,  and  the  error  in  the  approximation  \T{U)  — 

Figures  2 and  3 show  the  approximation  error  for  each  level. 

The  error  function  from  Level  1 clearly  demonstrates  the  ability  of  the  in- 
verse multiquadric  to  approximate  smooth  data.  The  error  near  the  boundary 
is  scaled  by  an  order  of  magnitude  at  each  level,  but  has  the  same  general  be- 
haviour. The  final  iteration  leaves  error  near  the  boundary,  aggravated  by  test 
points  being  outside  the  original  scattered  data  set.  Such  evaluation  points 
ought  to  be  included  since,  although  they  require  the  extrapolation  of  Sk  to 
evaluate,  the  experiment  was  specified  to  be  conducted  on  the  unit  square. 
The  original  aim  of  finding  a solution  where  the  residual  is  comparable  to  the 
approximation  accuracy  is  fulfilled  at  Level  3. 

Example  I is  repeated  as  far  as  the  regularization  of  data,  and  then  poly- 
nomial tensor  product  splines  are  used  to  find  the  final  approximation,  as 
described  in  Section  2.  Such  splines  cannot  replace  the  inverse  multiquadric 
approximation  on  the  regular  grid  without  an  increase  in  error.  Such  an  er- 
ror is  then  propagated  to  the  next  level  where  the  discrepancy  is  amplified. 
However,  if  the  hierarchical  refinement  procedure  is  abandoned,  then  these 
basis  functions  efficiently  yield  a good  approximation.  Table  2 shows  the 
approximation  accuracy  for  such  splines  of  different  polynomial  degree  with- 
out iterative  refinement.  The  grid  sizes  are  comparable  to  those  used  in  the 
original  example. 
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Grid  Size 

Linear 

Quadratic 

Cubic 

21  x 21 
41  x 41 
81  x 81 

2.2  x 1CT2 

5.3  x Hr3 

1.4  x nr3 

6.4  x nr3 
9 x 10~4 
1 x 10~4 

3.5  x HT3 

2.6  x HT4 

5.7  x 1(T5 

Tab.  2.  Approximation  error  for  the  various  splines. 


Fig.  2.  Approximation  Error  for  Levels  1 and  2. 


Fig.  3.  Approximation  Error  for  Level  3. 

§5.  Conclusion 

A global  solution  to  an  interpolation  problem  involving  a large  number  of 
data  points  is  too  expensive  to  compute  directly  if  inverse  multiquadrics  are 
to  be  used  effectively.  However,  if  the  aim  is  to  generate  approximations  to  a 
function,  then  such  a solution  is  unnecessary,  and  an  alternative  method  has 
been  presented. 

The  underlying  idea  is  to  transform  given  scattered  data  fi  at  points  Xi  to 
regular  approximate  data  fi  at  t/i,  which  is  easier  to  solve  for.  The  aim  is  then 
no  longer  to  interpolate  the  initial  data,  but  to  find  a good  approximation  to 
it.  The  final  solution  s is  an  approximation  to  /,  which  is  close  to  /.  Success 
relies  on  minimising  the  local  approximation  errors  | f{yi)  — f{yi)\- 

The  algorithm  is  0{N ) since  the  only  work  related  to  the  number  of  initial 
points  is  the  search  for  the  q closest  points  to  each  yi.  Such  a search  can  be 
improved  by  making  assumptions  as  to  which  Xj  are  unlikely  to  qualify. 
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The  time  required  to  solve  each  sub-domain  problem  is  reduced  due  to 
the  use  of  the  stored  n x n key  inverse  matrix.  Solving  directly  would  be 
0(n3),  but  the  required  matrix-vector  product  is  0(n2). 

The  hierarchical  iterative  refinement  strategy  produces  good  approxima- 
tions, and  is  the  only  sequential  aspect  of  the  method.  The  search  for  ap- 
proximate regular  data,  and  the  solutions  for  each  decomposed  sub-domain 
are  parallelizable  operations,  although  this  has  yet  to  be  implemented.  These 
features  mean  that  large  data  sets  can  be  dealt  with  in  acceptable  computing 
time. 
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